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Abstract 

Customer slowdown describes the phenomenon that a customer’s service requirement 
increases with experienced delay. In healthcare settings, there is substantial empirical 
evidence for slowdown, particularly when a patient’s delay exceeds a certain threshold. 

For such threshold slowdown situations, we design and analyze a many-server system that 
leads to a two-dimensional Markov process. Analysis of this system leads to insights into 
the potentially detrimental effects of slowdown, especially in heavy-traffic conditions. We 
quantify the consequences of underprovisioning due to neglecting slowdown, demonstrate 
the presence of a subtle bistable system behavior, and discuss in detail the snowball effect: 

A delayed customer has an increased service requirement, causing longer delays for other 
customers, who in turn due to slowdown might require longer service times. 

1 Introduction 

The phenomenon of customer slowdown describes the fact that a customer’s service requirement 
increases with the customer’s experienced delay. While the operations management literature 
is largely built on the assumption that service times are independent of delay, a growing number 
of empirical studies, predominantly in healthcare settings, provide evidence for situations where 
slowdown occurs. This empirical evidence calls for the development of stochastic models that 
take into account slowdown, in order to not only assess its impact on the performance of service 
operations, but also to gain understanding of the fundamental changes that slowdown brings 
to system behavior. 

A large body within the healthcare operations literature investigates the impact of workload 
on service times of patients. A canonical example in this domain is the admission of patients 
to the intensive care unit (ICU). There is substantial empirical evidence for slowdown in such 
settings: delays in receiving appropriate care can result in adverse effects such as an increased 
length of stay in the ICU 0 [HI |9l |22l |23l [H]. Since ICUs are typically heavily used and 
subject to unforeseen circumstances, delays in admitting patients are the rule rather than 
the exception, which makes the slowdown effect potentially threatening. A delayed patient 
that requires a longer service time will increase the overall workload of the system, therefore 
causing longer delays for other patients, who in turn due to slowdown might require longer 
service. This triggers a snowball effect, with an impact that is hard to assess without having 
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a detailed understanding of the global system behavior that takes into account the subtle 
dependencies among customers due to slowdown. Particularly when a system like an ICU 
is designed to operate under heavy-traffic conditions, the neglect of slowdown might lead to 
underprovisioning and severe performance degradation. 

Critical care systems such as an ICU are typically modeled as multi-server systems that 
operate in heavy-traffic regimes IllISl. The patients are the customers, the beds are the servers, 
and the performance analysis of the multi-server systems gives insight into the patient flow. 
We shall consider a Markovian multi-server system with the additional feature of slowdown. A 
detailed analysis of this system gives insight into the key features of slowdown, in particular 
when compared against multi-server systems without slowdown. 

1.1 A threshold slowdown system 

Slowdown can be modeled as a non-increasing function ^(•) that describes the rate of service 
as a function of the queue length seen upon arrival. That is, a customer meeting n customers 
upon arrival will receive service with rate /r(n), regardless of arrivals and departures after the 
customer has joined the system. Note that, since the number of customers seen on arrival can 
be translated into an expected delay, the service rate can also be interpreted as a non-increasing 
function of the expected delay. Assuming a service rate that is a function of the state of the 
system, leads to a so-called state-dependent queueing system. 

The majority of the empirical studies on slowdown has focused on a threshold slowdown: 
if a patient’s delay surpasses a certain threshold, he will receive a longer service time and 
otherwise he receives a service of regular length. In terms of the slowdown function, this means 
that fi{n) = fin n < N and fi{n) = if n > N with ^Uh > /^l- Note that the expected delay 
is translated to a number of customers n met on arrival and compared against the threshold 
N. The definition of the threshold varies across different medical conditions and situations. In 
[7] it is argued that a critically ill patient awaiting transfer from the emergency department to 
the ICU is labeled as delayed if the patient has waited longer than 6 hours. Delayed patients 
on average have an ICU length of stay that is 1 full day longer than the non-delayed average 
length of stay. Similar conclusions are drawn in [23] for the same situation in different hospitals. 
However, [23| uses a threshold of 8 hours. Both studies [Zl[23| establish a strong correlation 
between the delay a patient experiences in receiving an assigned bed and the ICU length of 
stay. Depending on the medical condition, the delay threshold can be in the order of minutes, 
such as for cardiac arrest patients |9|, hours, as seen in [3 EH], or even days, such as the 2 
day delay in receiving surgery [2Tj; or a 3 day threshold of delay for pneumonia patients |22j . 
An encompassing study is performed in |8|, where it is empirically verified that the slowdown 
effect is prevalent across multiple hospitals and patient conditions. 

We shall adopt the model in [8] , which is a multi-server model with a threshold service rate 
function ^(-j. Customers arrive according to a Poisson process with rate A, have an exponential 
service requirement, and are served by s servers. Due to the threshold, we then distinguish 
between two types of customers: those who were taken into service immediately upon arrival 
(non-delayed) and those who have experienced delay (delayed). We set the threshold to N = s 
so that non-delayed customers are served with a high service rate /in and delayed customers 
are served with a low service rate with //h > Indeed, in that case, delays cause a longer 
service time. Define the two-dimensional Markov process Y{t)), with X(t) G No the total 

number of customers in the system at time t, and Y(t) G {0,..., s} the number of non-delayed 
customers in service at time t. Denote = X/{spu) and Pl = A/(spl)- The system is stable 
when Pl < 1. When pn = Pl the model reduces to the standard M/M/s system. The load of 
the slowdown system is described by p = (1 — P(IU > 0 ))ph -|- IP(IU > 0)pl, where W is the 
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Figure 1: Transition rate diagram and state space of the threshold slowdown system. 


stationary waiting time, so that F{W > 0) is the delay probability. Figure [^displays the state 
space and the transition rate diagram. In |8], approximations are derived for key performance 
indicators that give insight into the slowdown effect. Based on parameter values calibrated 
from real ICU dataflows, the approximations in [8] indicate that the slowdown effect can be 
substantial, and should not be ignored in critical care systems that operate in heavy traffic. 

This two-dimensional Markov process can be used to investigate the impact of slowdown, 
both qualitatively and quantitatively, in particular in comparison with the widely applied 
M/M/s system (which neglects slowdown). While the focus in [S] lies on approximations for 
small to moderate-sized systems (ICUs of 6 and 15 beds), we focus on exact and asymptotic 
results, both for finite s and the regime s —>■ oo and p\ 1. With exact results we refer to deter¬ 
mining the stationary distribution of the Markov process using a numerically stable algorithm. 
This algorithm allows us to compute the exact two-dimensional stationary distribution for not 
only small but also large systems. Asymptotic results give rise to accurate approximations for 
the dimensioning of large systems in heavy traffic. Particularly, we are interested in the effect 
of slowdown in the Quality-and-EfRciency driven (QED) regime [16]. As it turns out, the way 
to establish non-degenerate limiting behavior for a multi-server system with slowdown in a 
QED-type regime is by letting Pl approach 1, and pn approach Pl as s —oo. We find that 
this scaling window is such that the probability of delay converges to a value that lies strictly 
in the interval (0,1), which is a manifestation of non-degenerate limiting behavior. As pointed 
out in |8|, deriving exact results becomes mathematically challenging because determining the 
stationary distribution of the Markov process involves high-dimensional matrix inversion. To 
relieve this computational burden of a large state space (particularly for large s), we exploit 
the fact that the Markov process has a block diagonal structure in the inner states (states with 
more than s customers in the system), which allows for an exact solution using matrix-analytic 
techniques. This technique typically relies on iterative algorithms that solve a non-linear ma¬ 
trix equation. For our model, we are able to find an exact solution for this matrix equation, 
which then immediately renders the problem of computing the stationary probabilities of the 
inner states computationally tractable, also for large s, see Section |4.1[ What remains is the 
computation of the stationary probabilities of the boundary states (states with s or less cus¬ 
tomers in the system). We introduce a novel approach that computes the exact stationary 
probabilities of the boundary states by exploiting the transition structure and by introducing 


first-passage probabilities. A detailed description of this approach can be found in Section 4.2 
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1.2 On the relation with operator slowdown 


Slowdown can refer to customer slowdown and operator slowdown. Customer slowdown refers 
to an increase of a customer’s service requirement, caused by the delay experienced by that 
customer. Operator slowdown refers to a service rate that decreases with the workload present 
in the system. Operator slowdown usually occurs in large service systems, such as call centers, 
due to fatigued operators [TO]. However, it is also common in medical applications under high 
workload, where care providers have to multitask and share (now crowded) central resources 
such as computer terminals |3|. The key difference is that customer slowdown starts with an 
individual delayed customer, and affects all customers behind this customer, while a decreased 
service rate in operator slowdown affects all customers that are in service. Customer slowdown 
therefore typically requires a more detailed state description, making it harder to analyze 
than operator slowdown. In this paper we indeed focus on customer slowdown, but we make 
comparisons with operator slowdown in several places. 

Operator slowdown under Markovian assumptions leads to a one-dimensional Markov pro¬ 
cess which is more tractable than our two-dimensional process and is amenable to fluid analysis. 
In |10] an M/M/s-type model with operator slowdown is investigated. Additional properties 
in m are customer abandonments and state-dependent service rates. We make a comparison 
with [TO| by extending our base model to also include customer abandonments in Section 2.2 


Both slowdown models exhibit a bistable behavior in which the models alternate between two 
dominant regions. For the customer slowdown model, however, this behavior is more subtle 
than for the operator slowdown model (see Section 2.2). 

Both customer and operator slowdown fall into the broad category of queueing systems with 
state-dependent service rates, like for instance an MjGjl system with state-dependent service 
rates m- In 1^, the optimal admission policy is studied for an MjGjX system with service 
rates that increase with the workload below a certain threshold and decrease with the workload 
above this threshold. State-dependent queueing systems also arise when arrival and/or service 
rates are dynamically controlled to minimize average cost per time unit, see e.g. [SlITOlllHI. 
All these examples concern operator slowdown. 


1.3 Structure of the paper 

The paper is structured as follows. Based on a detailed analysis of the two-dimensional Markov 
process in Figure we identify three key features of threshold slowdown systems; severe 
performance degradation due to the snowball effect; a subtle bistable system behavior; and the 
existence of non-degenerate limiting behavior in a QED-type heavy-traffic regime. We discuss 
these three features in Section]^ The first two features were identified by using the stationary 
distribution of the two-dimensional Markov process. Section introduces the model in greater 
detail and Section describes how we solve for its stationary distribution using matrix-analytic 
methods and some properties of regenerative processes. The QED-type heavy-traffic regime is 
outlined in Section We conclude in Section and present some supporting results in the 
appendix. 


2 Key features of threshold slowdown systems 

Unless stated otherwise, we assume a stable system, i.e. pi^ <1. 
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Figure 2: Sample path of the total number of customers in the system L{t) as a function of 
time t, for s = 15, A = s and loads pu = 0.7 and Pl = 0.98. 

2.1 Performance degradation due to the snowball effect 

We first present a detailed description of the snowball effect caused by slowdown and then 
assess the adverse effects for system performance. 

For explanation purposes, we refer with busy periods and idle periods to the excursions 
of the process X(-) above and at level s, and below level s, respectively. Hence, during busy 
periods, newly arriving customers will experience delay and are thus subject to slowdown. The 
snowball effect sets in each time a new busy period starts. An example sample path of idle and 
busy periods is given in Figure where we plot the total number of customers in the system 
L(t) at time t. Compared with an M/M/s system without slowdown (with a high service rate 
Pu), the busy period in the time interval (9300,10000) is relatively long, due to the slowdown of 
delayed customers that reinforces, through other delayed customers, the persistence of the busy 
period. Such busy periods are essentially equivalent to busy periods in an M/M/s system with 
a low service rate /Tl- These excursions during which congestion levels are high occur relatively 
frequently due to the snowball effect that triggers them, and this leads to severe performance 
degradation, particularly in heavy traffic. 

This performance degradation is visible in Figure]^ which displays for the same parameter 
values as in Figure the stationary distribution of the total number of customers in the 
system L of the threshold slowdown system. This stationary distribution is calculated using the 
numerical scheme that will be discussed in Section We also plot the stationary distribution 
of an M/M/s system with uniform service rate pn (the fast system) and with uniform service 
rate pi^ (the slow system). We append the subscripts H, or L to random variables to indicate 
that they belong to the M/M/s system with high service rate, or low service rate, respectively. 
We see that the distribution of the slowdown system peaks around the same point as the fast 
system, but that the tail behavior of the slowdown system is more comparable to the slow 
system (which can be attributed to the snowball effect and long busy periods). Such a fat tail 
obviously has severe consequences for performance, and in Figure we see for instance that 
the mean number of customers in the system increases considerably due to slowdown. 

The neglect of slowdown might lead to underprovisioning. Table provides an example 
in which we search for the number of servers s that are required to achieve a certain delay 
probability. Naturally, the threshold slowdown model requires equally many or more servers 
as required by the fast system with uniform service rate ^h- In particular, differences between 
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Figure 3: Stationary distributions of the total number of customers in the system L for systems 
with and without slowdown. Parameter values are s = 15, A = s, /Oh = 0.7 and Pl = 0.98. The 
expected number of customers in the systems are IE[Lh] ~ 10.8, IE[Ll] ~ 59.4 and IE[L] ~ 34.5. 


Mh pL A 

F{W > 0) 4 s* 

F{W > 0) si s* 

1 0.9 10 

0.1 16 16 

0.5 12 13 

12 

18 18 

14 15 

15 

22 22 

18 19 

20 

27 28 

23 24 

1 0.7 10 

0.1 16 17 

0.5 12 15 

12 

18 19 

14 18 

15 

22 23 

18 22 

20 

27 30 

23 29 


Table 1: Minimal number of servers (fast system) and s* (threshold slowdown system) 
required to achieve a certain P(VF > 0). 

the required number of servers in the fast model and the threshold slowdown model seem to 
increase with the delay probability, with the ratio /Uh///l and with the arrival rate A. 

Another indicator for substantial slowdown effect is the difference p — pn, as we will show 
in the next subsection. This difference is the increase in load caused by the slowdown effect 
with respect to the load of the fast system. 

2.2 A subtle bistable behavior 

The threshold slowdown system behaves as the fast system below the threshold, and as the 
slow system above the threshold. However, for many relevant parameter settings, neither the 
fast nor the slow system provides a good approximation for the slowdown system. The reason 
is that the slowdown system in fact is a rather intricate mixture of both system as will be 
explained in this subsection. 

We start by examining the two-dimensional stationary distribution, which typically consists 
of two dominant regions; region 1 with only non-delayed customers and no customers waiting, 
and region 2 with delayed (slowdown) customers in service only and many waiting customers. 
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(a) Ph = 0.6, p- pn^ 0.123 and ¥{W > 0) = 0.32. 



s 2s 3s 

(b) Ph = 0.95, P - Ph = 0.027 and P(W > 0) = 0.90. 


Figure 4: For parameter settings that are mild or extreme the slowdown system resembles 
either the fast or the slow system. The horizontal axis plots the total number of customers in 
the system and the vertical axis indicates the number of non-delayed customers in service. The 
contour plot shows where the probability mass is located (darker colour means more mass). 
Parameter values are s = 15 and Pl = 0.98. 


Region 1 thus complies with the fast system and region 2 with the slow system. An important 
parameter that determines whether region 1 or region 2 is dominant is pn- A low to moderate 
Ph makes region 1 dominant, which suggests using the fast system as a proxy. A high load Ph 
makes region 2 more important, and in fact, when pn approach 1, the slow system will be a 
good approximation. See the two examples in Figure Notice here that for a system with a 
high delay probability, i.e. Figure Qb), the increase in load p — Ph due to the slowdown effect 
is small, since both loads pn and Pl are large and comparable. In contrast, the increase in load 
in Figure Qa) is much larger. 

Arguably the most natural scenario, when pn is high but not extremely high, say pn € 
(0.7, 0.9), gives a less clean picture. Then the slowdown system is a mixture of the fast and 
slow systems, under the right condition that Pl is decisively larger than pn. A good example 
is Ph = 0.8 and Pl = 0.98, as can be seen in Figure [^b). By increasing the load Pl busy 
periods become longer, causing the shift in probability mass towards region 2 and increasing 
the severity of the slowdown effect in terms of p — pn as is witnessed in Figures [^a)-(b). 

Our slowdown system thus has a subtle bistable behavior, which rises to the surface when 
both Ph and Pl — Ph are substantial but not extreme. A more extreme bistability effect would 
occur when pn could become larger than 1. We therefore next discuss two extensions of our 
model that allow for Pl > 1: 

(i) A threshold slowdown system with a finite waiting room; 

(ii) A threshold slowdown system with customer abandonments. 

System (i) can have at most N customers in the system and is therefore inherently stable. 
When Pl > 1 and pn is sufficiently small the system will alternate between periods during which 
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(a) Pl = 0.9, p - pH = 0.047 and P(H/ > 0) = 0.46. 



s 2s 3s 

(b) Pl = 0.98, p- pn= 0.144 and ¥{W > 0) = 0.80. 


Figure 5: Two dominant regions in the stationary distribution become visible when the load 
of the delayed customers /?l increases. Parameter values are s = 15 and pu = 0.8. 


the process settles in the high-occupancy states around N, and periods in low-occupancy states 
below s. This gives rise to bistable behavior, and for some parameter ranges even leads to a 
bimodal distribution as seen in Figure j^a). This bimodality can be explained by the fact that 
for Pl > 1; the process has two clear points of attraction: the state N and the state p^s where 
the rate of arriving and departing customers is equal. Note that our original slowdown system 
has only one point of attraction, because pu < pi^ < 1. 

System (ii) assumes that waiting customers abandon the system after an exponential time 
with mean 1 /6. Because the total abandonment rate scales linearly with the number of waiting 
customers, also this system is inherently stable. For Pl > 1 it has two points of attraction: one 
below s, and one above s precisely where the total rate of arriving customers equals the rate 
of departing (abandoning and served) customers. For Pl > 1 this process alternates between 
the two points of attraction as is shown in Figure j^b). This system is closely related to the 
operator slowdown system with abandoning customers considered in [lUj . In jlU| , the bistability 
effect was also observed, where the two points of attraction were identified explicitly. Explicitly 
characterizing the two points of attraction in the customer slowdown model is more difficult 
due to the two-dimensional nature of the system. 

2.3 Scaling limits 

So far we increase either the number of servers or the arrival rate. We continue by examining 
a scaling of both parameters at the same time. It is well known in the literature that for 
GjMjs queues, one should scale the arrival rate or the number of servers such that the load 
p(^) ~ 1 — [3/^/s as s —7- oo to achieve QED performance [l6]. In terms of our model parameters, 
the scaling is then as follows: 

(3/^/s), ( 2 . 1 ) 
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(a) A finite buffer that limits the total number (b) Allowing customers to abandon the queue 
of customers in the system to N. Here, s = 81, with rate 5. Here, s = 36, pn = 0.7 and pL = 
N = 93 and pH = 0.8. 1.2. 


Figure 6 : Two extensions of the multi-server system with slowdown that exhibit, for a narrow 
range of parameter values, a bistability effect that is visible in the bimodal marginal distribution 
of the total number of customers in the system. 


with constant (5 > 0 and s > to guarantee a positive arrival rate By applying (2.1) 


to our multi-server system with slowdown one finds that we establish so-called Quality-Driven 
(QD) performance. QD performance refers to a very high quality of service, e.g. the probability 
of delay goes to 0 and many servers are idle. This might be undesirable in view of unnecessary 
operational costs (overdimensioning). The reason for QD performance is that since /Th > /^l 
we have < 1 in the limit for s —)• oo. This ensures that the system stabilizes around a 
state with relatively low occupancy and with only non-delayed customers in service and no 
customers waiting in the system. To obtain QED system behavior we set the high service rate 
according to 


= fiL(l +7/Vs), 


( 2 . 2 ) 


with constant 7 > 0. Note that now /Xh^//xl —)• 1 for s —f 00 and thus also goes to 1. We 


refer to the combination of ( 2 . 1 ) and 


as a QED-type regime. The reason for this choice 


of scaling becomes clear when we examine the load of the slowdown system with s servers 


- 


= 1 -^ 


p . 1 + P(W(") > 0)^ 




1 + ^ 


(2.3) 


which shows that t 1 s-s s —f 00 . Compared to the standard scaling of the load in GjMjs 
queues, the load in the customer slowdown model approaches 1 slower as it is multiplied by 
the second term in (2.3). Figure depicts the probability of delay as a function of s, which 
indeed shows that the probability of delay converges to a value in ( 0 , 1 ). 

Using stochastic coupling techniques, we related the two-dimensional process to two one¬ 
dimensional processes that serve as a stochastic lower and upper bound (at the process level 
in terms of stochastic domination). These two related processes are the fast and the slow 
system introduced earlier. The bounding processes provide sharp approximations for the two- 
dimensional process. For both bounding processes, we show that in the QED-type regime, the 
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Figure 7: Probability to wait for the fast, slow, and customer slowdown system. For all three 
systems, the scaling (2.1) and (2.2) is used with (/3,7) = (0.5, 0.5). 


probability of delay converges to a value strictly in between 0 and 1 , and this then also holds 
for the two-dimensional process. Hence, this provides strong evidence for the existence of a 
non-trivial stochastic-process limit. Formally establishing the existence and characterizing this 
stochastic-process limit is a challenging open problem, because the limiting process is likely to 
be two-dimensional as well, and classical techniques to prove stochastic-process limits [ 2 U] do 
not seem to carry over easily. 

2.4 Insights 

Here we summarize the insights obtained in this section. 

Customer slowdown of the threshold type leads to severe performance degradation, particu¬ 
larly in heavy traffic. Compared to a system without slowdown, the busy periods are relatively 
long due to the slowdown of delayed customers that reinforces the persistence of the busy pe¬ 
riod. We refer to this effect as the snowball effect, which describes the correlated service times 
when customers are delayed. Further, for a relatively high load pn, we find that the threshold 
slowdown system is a mixture of the fast and slow systems. This mixture effect is visible in 
the two-dimensional stationary distribution, where it manifests itself as two dominant regions 
in terms of probability mass - a subtle bistable effect. Finally, by using a QED-type scaling 
for the arrival rate and the fast service rate pn; we have shown that a non-degenerate limit 
behavior occurs as the number of servers increases. 

3 Model description 

Recall that X{t) G No is the total number of customers in the system at time t and Y{t) G 
{0,1,..., s} is the number of non-delayed customers in the system at time t. Note that X{t) > 
Y(t). Then, {(X(t), Y{t)), t > 0} is an irreducible continuous-time Markov chain with discrete 
state space V UW, given by 

V = {{hj) \0 <i < s, 0 < j <i}, W = {{i,j) \ i > s, 0 < j < s}. (3.1) 
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Recall that we refer to the states with s or less customers in the system, as the boundary states. 
With inner states we refer to the states with more than s customers in the system. For an 
inner state (i,j) with i > s, we have the following transition rates: 

• From (i, j) to (i + l,j) with rate A, 0 < j < s; 

• From (i,j) to (i — l,j) with rate (s — 0 < j < s; 

• From (i,j) to (i — l,j — 1) with rate I < j < s. 

The transition rate diagram of the continuous-time Markov chain is shown in Figure 

Define level i as the set of all states with a total of i customers in the system. Now we 
have the following alternative description of the transition rates. The matrices contain the 
transition rates from level i to level i + k with i > s. Let I be the identity matrix of size s -|- 1. 
Then the matrices are given by Ai = A/, 


Ao 


f X + sni^ 

A -|- (s — 1 )//l + /^H 


\ 


V 


A + S/in/ 


(3.2) 


and 


/sfJ-L 

Hh {s - 1)^L 


A-i = 


2/Xh 


V 


\ 


SHu 0/ 


(3.3) 


By assumption < pi^, and we have the following condition for ergodicity of the Markov 
process. 

Lemma 3.1. The Markov process is ergodic if and only if 


Pc < 1. (3.4) 

Proof. We require that the mean drift in the negative direction is larger than the mean drift 
in the positive direction; see Neuts’ mean drift condition 1211 Theorem 1.7.1]. This gives 

ttAiI < 7rA_il, (3.5) 

where 1 is a column vector of ones of size s -|- 1, tt is the solution of tt X)fc=-i Afc = 0 with 
ttI = 1. We clearly have tt = (1, 0,..., 0) and thus the result follows. □ 

4 Obtaining the stationary distribution 

Assume that ( |3.4[ ) holds and define the stationary probabilities 

p{i,j):= \imF{X{t) = i,Y{t)=j), {i,j)eVUW. (4.1) 
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The balance equations for the inner states are obtained by equating the rate out of and into 
an inner state (i, j), yielding, for z > s, 0 < j < s, 


(A ++ {s - j)l^i.)p{i,j) = Xp{i - l,j) + {s-j)pLP{i + l,j) 

+ (j + 1 )^hP(* + l,j + 1), 


(4.2) 


where by convention p{i,s + 1) = 0. Equations (4.2) are referred to as the inner equations. 


The balance equations for states with i < s are called the boundary equations. 

The stationary probabilities of the inner states are determined using matrix-analytic meth¬ 
ods that search for the solution to a non-linear matrix equation. Exploiting structural proper¬ 
ties of the Markov process, we derive explicit solutions for these matrices. Eor similar explicit 
results using the matrix-analytic methods, see 1251 [261 127]. Next, we solve the boundary equa¬ 
tions. Since we want to be able to solve the stationary distribution also for large s, solving the 
(s-|-l)(s-|-2)/2 boundary equations using Gaussian elimination might become computationally 
cumbersome. We therefore present a more sophisticated algorithm that exploits the structure 
of the state space and the explicit matrix solution. 


4.1 Inner equations 

Let Pi = (p(z, 0),p{i, 1),... ,p{i, s)), and rewrite the inner balance equations as 

Pi_iAi -h PiAo + Pi-HiA_i = 0, i> s. (4.3) 

The rate matrix R is defined as the minimal non-negative solution of the non-linear matrix 
equation (H] Theorem 3.1.1] 

Ai + RAo + R^A_i = 0. (4.4) 

It can be shown that the stationary probabilities satisfy 

Pi = PsA’*“®, i > s. (4.5) 


Since the transition matrices are all lower triangular, so is the rate matrix R. Denote 

/ -^0,0 \ 

Ri,o Ri,i 


R = 


\^s,C 


(4.6) 


Rs,s/ 


and note that R^ is again a lower triangular matrix with elements {R^)ij = Yl\=j Ri,kRk,j for 

* > j- 

Equation (4.4) consists of (s -(- 1)^ separate equations. Eor the diagonal elements we have 


A - (A -h (s - j)/ZL + jPH)Rj,j + {s- j)Pi.Rj,j = 0, 0 < j < s, 

A — (A -|- spii)Rs^s = 0, j = s, 


(4.7) 

(4.8) 


where Rjj in (4.7) is obtained as the minimal non-negative solution. The minimal non-negative 
solution of (4.7) is contained in the interval (0,1), because for Rjj = 0 the left-hand side of 
(4.7) is positive, for Rjj = 1 the left-hand side of (4.7) is negative, and we are dealing with 
a continuous function. Interestingly, Rq q = Pl and Rjj is monotonically decreasing in j. Eor 
each element on the subdiagonals we have a linear equation with solution 


Ri,j ~ 


X)fc=j-|-1 J (s + Y2k=j+1 ^iA^kJ+liJ A pH 

-^ + ('S — j)pL + jpn — {Ri,i + Rj,j) (s ~ j)/^L 


(4.9) 
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foT j = i — h, h < i < s and h = 1,2,..., s. In (4.9) we use the convention that /(^) ~ 


if ii < io- Equations (4.7)-(4.9) fully describe the rate matrix R. 


Recall that a lower triangular matrix is non-singular if it has all non-zero elements on the 
diagonal. Thus, the matrix R is non-singular and also I — R is non-singular. The inverse of 
I — Ris required to compute the stationary probabilities, as the normalization of the stationary 
distribution partially follows from p 5 (/ -|- i? -|- -|- • • • )1 = Ps(/ — R)~^l. The elements of the 

inverse are given by 


((/ - R)-%j = 


{{I-R)-%j = 


1 


(I-R) 


0 < j < s. 




{I - R)i,i 


, 0 < j < i < s. 


(4.10) 

(4.11) 


Instead of searching for R, one can also search for the matrix G, defined as the minimal 
non-negative solution of the non-linear matrix equation 


A_i + AoG +AiG^ = 0. (4.12) 

The matrices R and G are related as AiG = RK-i and thus G = A)”^RA_i, which exists since 
Ai is a diagonal matrix. 


4.2 Boundary equations 

The boundary equations can be represented as a set of (s -|- l)(s -|- 2)/2 linear equations, which 
can be solved using Gaussian elimination with an arithmetic complexity of 0(s®) |12l Chap¬ 
ter 10]. By exploiting the structure of the boundary equations one can reduce the arithmetic 
complexity to O(s^). In short, we wish to embed the Markov process on level s for which we 
need the return probabilities when jumping to a higher level (the matrix G), combined with 
the return probabilities when jumping to a level below (yet to be determined). This allows 
us to compute the non-normalized stationary probabilities of the states in level s. Then, we 
recursively compute the remaining boundary probabilities starting from level s —1, working our 
way down to level 0, leading to a non-normalized stationary distribution. Finally, the normal¬ 
ized stationary distribution follows by summing over all states and dividing the non-normalized 
version of the stationary distribution by the resulting sum. 

To this end we introduce two first passage probabilities. To aid the derivation of these 
probabilities we introduce subsets of V, indexed by a state (k, 1) G V. We identify the triangular 
set of states T(^k,i) fo the South-West of the state {k, 1), specifically, T(^k,i) •= {(b j) \ k — I < i < 
k — I, 0<j<i — {k — Z)}, see Figure]^ 

Let 0k(i,j) be the first passage probability to phase j -|- 1 in state {i + 1 — k, j + 1), where 
the Markov process starts in state (i,j) G T(^s,s-k)- Note that by phase j we refer to the set of 
states with Y(t) = j. By one-step analysis we obtain 


^o(bj) 
6'fc(b j) 


A • 1 -|- (i 


A • 0 -h (f 


j)ML • 0 + j>H6'o(i - l,j - l)<9o(i,j) 

A+ (* - j)l^L + JAiH 

j)^L6*fc-i(i - 1, j) +J>HEf=o - 1 

A + (« -i)AtL +JAh 


(b j ) ^ ) 

j - l)0k-i{i - l,j) 


(bj) ^ 'R'{s,s—k)'! Zc > 0, 


(4.13) 


(4.14) 
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Yit) 



{k,l) 


Xit) 


Figure 8: A visual clarification of the triangular set of states 


from which the following expressions can be readily derived, 


^o(bj) r "i /-■ -V "i' ■ /I 1 ' 1 ’ (bj) ^ 

A + (* - + J/^h( 1 - 6'o(* - 1,J - 1)) ^ ^ 

Of _ (i - j)k-Ldk-i{i - 'i-J) + - l,j - l)0k-iii - l,j) 

“ A + (i - j>, + jMH(l - 0o(i -l,j- 1)) 

(bj) e 7^(5,s-fc), > 0. 


(4.15) 


(4.16) 


Note that 9o{i,i) = 1, which means that if the Markov process reaches a state on the main 
diagonal, it eventually always reaches state (s, s). 

Let ip[k,i){hj) be the first passage probability to level k in state (fc,/), where the Markov 
process starts from state (i,j) G We express these hrst passage probabilities in terms of 

^kihj) follows 


i+l 

i’{k,i){hj) = + {i,j) eT(^k,i)- 


(4.17) 


The computation of the hrst passage probabilities ^|J is the most time consuming step in the 


derivation of the boundary probabilities. Equation (4.17) is evaluated for a total of s(s + l)(s + 
2)(s + 3)/24 combinations of {i,j) and {k,l), leading to an arithmetic complexity of O(s^). 

Let 'L be the matrix of elements '4’{s,k)i^~ 1) j) where j is the row number and k the column 
number. The balance equations of the Markov process embedded at level s are 


i(A_i'I' + Ao + AiG) = 0. 


(4.18) 


One solves this homogeneous set of equations using the numerically stable Grassmann version 
of the Gaussian elimination algorithm [TT] to obtain the stationary probabilities at level s. This 
solution is not normalized. 

In order to obtain the remaining boundary probabilities, one embeds the Markov process 
on the levels i, f + 1,..., s,... with i < s, for which we derive the following balance equations: 

+ (* - j)/^L + jMH(l - 0o{i -l,j - 1))) 

= P{i + 1, j)(* + 1 - J>L +p(* + 1, j + l)(j + 1 )Mh 

+ '^P{hk){i - - l,k) + ^p(i,/c)fe/iHV’(i,j)(* -l,k- 1), j < i. (4.19) 

fc=0 k=l 
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One recursively solves (4.19) by first computing p{i,0), followed by p{i,l), et cetera, until 


p{i,i) is computed. Then, the remaining boundary probabilities follow by solving (4.19) for 
i = s — l,s — 2,...,1. Finally, the probability of having an empty system is found by examining 
the balance equation in state (0, 0), so that 


p(0,0) = (p(l,0)^L +p{l, 1)Mh)/A. 


(4.20) 


Recall that the stationary probabilities of level s are not normalized. Thus, the obtained 


stationary distribution of the boundary and inner states (obtained through (4.5)) are yet to be 


normalized. The normalized stationary distribution follows by dividing each non-normalized 
stationary probability by the sum over all states p(*, j) + Ps(.?^ — R)”^l. 

Using the stationary distribution, one can now obtain performance measures such as the 
delay probability 


P(1U > 0) = Ps+jl = Ps{I - R) ^1, 


(4.21) 


1=0 


or the mean queue length 


= ^*Ps+il = PsR(-f - R) ^1- 


(4.22) 


i=0 


4.3 Extensions 

We next describe how to obtain the stationary distribution of the slowdown model with (i) a 
finite buffer; or (ii) customer abandonments. 


4.3.1 Finite buffer 

The transition rate diagram of the slowdown system with a finite buffer is identical to the one 
shown in Figure but truncated at level N. Recall that we defined the matrices to contain 
the transition rates from level i to level i + k; since we now introduced the finite buffer, we 
restrict i as s < i < N. At level N, the matrix containing the transitions rates to level N — 1 
remains unchanged and is still A_i. The only difference is that there are no jumps in the 
positive direction and thus Ai does not exist and therefore the matrix Aq changes at level N, 
now indexed by an additional subscript N and given by Aq^at = Aq + Ai. 

The equilibrium equations in vector-matrix form are given by 

Pi_iAi PjAo + Pi+iA_i = 0, s <i < N, (4.23) 

Ptv-iAi + pArAo,iv = 0. (4.24) 

We now have the following relation, see m Section 2.2], 

Pi = pi-iRi, s <i < N, (4.25) 


where Ri is a level-dependent rate matrix with identical interpretation as the standard rate 
matrix of the matrix-geometric approach. One can now solve for the rate matrix by 


manipulating (4.24) as follows 


PiV = —PAf-lAi(Ao,Ar) ^ = pn-iRn- 


(4.26) 
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Note that Aq^tv is a diagonal matrix with negative elements and is therefore indeed non-singular. 
The remaining rate matrices are found from (4.23) as 


Pi =-pj_iAi Ao + i2j+iA-i =Y>i-iRi, s<i<N. 


-1 


(4.27) 


The matrix Aq + i?i+iA_i is lower triangular with negative elements on the diagonal and is 
therefore non-singular; for the proof of this statement, see [5l p. 519]. 

This leaves us to compute p^ and the equilibrium probabilities of the boundary states. We 
do so with the approach we have derived earlier for the slowdown system with an infinite buffer. 
The missing ingredients are the first passage probabilities from level s -|- 1 to level s, which are 
found through the relation 

Gi = A]'ii?iA_i, i > s. (4.28) 


Note that the auxiliary matrices Gi are level-dependent and have the same interpretation 
as the standard auxiliary matrix in the matrix-analytic approach. Thus, we substitute the 
level-dependent matrix G^+i for G in (4.18) and are able to compute the complete stationary 
distribution. 


4.3.2 Customer abandonments 


The base model is appended by adding transitions with rate 16 from state (s + I, j) to state 
{s + I — l,j) for I > 0. These transitions model a waiting customer abandoning the queue. By 
appending the base model with these transitions a level-dependent QBD (LDQBD) process is 
created. We use solution techniques for LDQBD processes as presented in UM to compute 
the stationary distribution of the semi-infinite strip of states and once again use the earlier 
derived technique to compute the equilibrium distribution of the boundary states. We briefly 
sketch the solution approach here. 

The aggregated abandonment rate depends on the number of customers waiting in the 
queue. This leads to level-dependent transition rate matrices which we label with an additional 
subscript I, such that describes the transition rates from level s -|- / to level s + I — k with 
I > 0. The transition rate matrices are given by Ai^; = Ai, Ag,; = Ao—161 and A_i^; = A-i+161. 

The solution approach is based on the same premise as for the finite QBD process case, 
namely 

Pi = Pi-iRi, i > s, (4.29) 


where Ri is a level-dependent rate matrix with identical interpretation as the standard rate 
matrix of the matrix-geometric approach. 

The following is explained in greater detail in [5]. Since generally only numerical solutions 
can be found for the Ri matrices of LDQBD processes, one resorts to computing the sequence 
{Ri}s<i<N*, where N* is chosen “large enough”. By [Sj Lemma 1] we have the explicit expres¬ 
sion 

oo j-i 

= E ^ n (4.30) 

j=0 k=0 


where Uf and are matrices defined recursively and are a function of the level-dependent 
transition matrices. Truncating the infinite sum in (4.30), one computes Rj^*. The remaining 
rate matrices then follow from the standard relation 


Ri — Ai (Ag^i—s “1“ A_i 


-1 


^ > s. 


(4.31) 


Note that the inverse exists. 
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As in the finite buffer case, this leaves us to compute ps and the equilibrium probabilities 
of the boundary states. Once again, the hrst passage probabilities are needed and follow 


from (4.28). Then, we substitute the level-dependent matrix for G in (4.18) and are able 


to compute the complete stationary distribution. 


5 A QED-type regime 

We next analyze the behavior of the multi-server queueing system incorporating slowdown for 
large s and Pl —)■ 1 by considering a sequence of queues, indexed by s. We write {X{t),Y(t))) = 
A = /Th = and pi^ = Without loss of generality we keep /Tl 

is) 

constant and assume throughout that > A^l- 

Let > 0) denote the probability that a customer has to wait in a slowdown system 

with s servers. We will identify a regime in which > 0) —)■ P(W > 0) G (0,1) so that 

the limiting system displays non-degenerate behavior, as in the classical QED regime. In order 
to do so, we introduce a random variable X^^^ (t) that represents the total number of customers 
at time t in an M/M/s queue where all customers are served with the high service rate /Uh^. 
As we have done before, we refer to this queueing system as the fast system. Let the random 
variable x[^\t) represent the total number of customers at time t in an M/M/s queue where 
all customers are served with the low service rate /Tl- We refer to this queueing system as the 
slow system. 

For two real-valued random variables A and B, we say that A stochastically dominates B 
if 

P(A < x) < F{B < x), (5.1) 

and we denote this as AY B. The following result is proved in Appendix [X| 


Lemma 5.1. x[^\t) Y Y xj/^(t). 

We next introduce the scaling 

A^®) = s/Xl - /3/XlVs, 
+ 7/Vs), 

with constants /Xl,/ 3,7 > 0 and s > /3^. Note that 




(5.2) 

(5.3) 

1 for s —)> oo. We refer to 
the scaling (5.2) and (5.3) as a QED-type scaling regime. We introduce the scaled random 
variables 








5.1 


Note that Lemma 

(t). The following lemma is proved in Appendix [B 


also holds for these scaled ranc 


oW(*) := 

Vs 

is). 


(5.4) 


om variables, i.e. dV [t) Y Y 


Lemma 5.2. If = dy'’ and Dh( 0) = a.s. with d]f’ —)■ d^, and Dl'’''(0) = df/’ and 

DVV — V a.s. with d\/'^ —)■ di^, then for s —)■ oo, and for every t > 0, the scaled random 
variables Dlf\t) V D^{t) and D[f\t) V Di^{t), where the infinitesimal means of the diffusion 
processes are given by 

a: < 0, 

^^l(-/3-7), X > 0, 

and constant infinitesimal variances (t|(x) = crL(x) = 2p,i^. 




(s) 




M) 


m^{x) = 


mi^{x) = 


Pi^{-l5-x), X < 0, 
-/3/Xl, X > 0, 


(5.5) 
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Remark 5.3. Both processes D^{-) and behave as an Ornstein-Uhlenbeck process below 
level zero and a reflected Brownian motion above level zero. 

Corollary 5.4. The probability density functions of Dii{oo) and Dl(oo) are given by 


fonix) 


0 ( 3 :+/ 3 + 7 ) 
'--H $(/3+7) > 




X < 0, 
X > 0, 


fD^{x) 


, X ^ u, 

(l-C0/3e-/^^ x>0, 


(5.6) 


with 


C„ = 


/3 + 7 


13 + 7 + 


4>{0++ 

$(/3+7) 


a. = 


/3 


f3 + 


m 

$(/3) 


(5.7) 


and (j){x) and <l>(x) the probability density function and cumulative density function of the 
standard normal distribution. 


Proof. Since we are dealing with piecewise-linear diffusion processes, one computes the proba¬ 
bility density functions using [6l Sections 1 and 4]. □ 


Remark 5.5. The stationary distribution of the diffusion process related to the fast system 
is equal to the distribution of the limiting random variable of the sequence (oo) — s)/\/s 
as shown in m Corollary 2], which establishes that an interchange of limits is allowed. Thus, 
one can use Xh'^^(oo) « s -|- Dii{oo)^/s. The same applies for the slow system. 

Corollary 5.6. The limiting probability of delay in the slowdown system P(1T(^) > 0) —)• 
P(1T > 0) G (0,1) for s oo and can be bounded as follows 

^ ^ > 0) = (l + (/? + 7)||^^)"'. (5.8) 


Proof. The limiting probability of delay in the fast system is computed from the distribution 
of Dii{oo) as 

roo 

P(1Th > 0) = / fouix) dx = 1 - Ch G (0,1) (5.9) 


and identically for the slow system to get P(1 Tl > 0) = 1 — Cl G (0,1). Using Lemma 5.1 


we find that these are lower and upper bounds on the limiting probability to wait ¥’{W > 0), 
respectively. □ 


6 Conclusion 

We have studied a threshold slowdown system in a Markovian setting. The threshold slowdown 
system incorporates a slowdown effect in which customers that are delayed require a longer 
service time. A delayed customer requiring a longer service time will increase the overall 
workload in the system, therefore causing longer delays for other customers, who in turn due to 
slowdown also require a longer service time. We refer to this phenomenon as the snowball effect. 
The snowball effect has been shown to be the leading cause of a severe performance degradation 
and the neglect of slowdown might lead to underprovisioning. A subtle bistable behavior is 
witnessed for slowdown systems with relevant parameter settings: the slowdown system either 
has only non-delayed customers and no customers waiting, or only delayed customers with 
many customers waiting, switching between configurations over time. We have introduced a 
QED-type regime for the slowdown system with many-servers and established non-degenerate 
limiting behavior. The snowball effect has been shown to persist in this QED-type regime. 
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A 


Proof of Lemma 


5.1 


The proof is based on a coupling argument and follows the same reasoning as [H Appendix B]. 
We distinguish between two cases: (i) ^ and (ii) X^^\t) ^ X^\t). Recall 

that for two real-valued random variables A and B, we say that A first-order stochastically 
dominates B if 

P(A < x) < P(i? < x). (A.l) 

(i) Assume that both queues see a common arrival process. Let the service time for the f-th 
arriving customer in the slow system be Bi^[i)] the corresponding service time in the slowdown 
model is then either B{i) = Bi^{i) or B{i) = Bi^{i) depending on whether the slowdown 

model has high (s or more customers in the system) or low congestion (less than s customers 
in the system) upon arrival of the f-th customer. Finally, we assume that both systems start 
empty. Let W{i) and Wi^{i) denote the waiting time of the f-th arriving customer before 
beginning service in the slowdown and slow system, respectively. We have the following result. 

Lemma A.l. Wl (i) > W{i) a.s. for all i. Moreover, xlf\t) > X^^\t) a.s. for all t. 


Proof. Let us prove the first statement using induction and fix an arbitrary event in the sample 
space a; G O leading to a sample path of the process. We append the argument u to the 
variables to indicate that we are studying a fixed sample path. Since we start with an empty 
system, observe that for the first customer we have = W{1,uj) = 0. Assume that the 

statement is true for the j-th arriving customer and consider the (j -|- l)-th arriving customer. 
For the sake of contradiction assume Wi,{j -|- l,w) < W{j + When customer j + 1 starts 

service in the slow system: 

• There are at most s — 1 customers among the first j arriving customers present in the 
slow system. 


• At least s customers from among the first j arriving customers are still present in the 
slowdown system since customer j + l has not yet started service in the slowdown system. 


From these two observations we conclude that there is a customer among the first j arriving 
customers that finished service strictly earlier in the slow system than in the slowdown system. 
However, due to the coupling we have B{i,uj) < Bi^{i,u}), i = and thus we have a 

contradiction. We have consequently established that Wi^{i,co) > W{i,uj) for all i. Recall 
that we fixed an arbitrary event and thus it holds for all cj G H. The latter statement of the 
proposition follows immediately. □ 


Lemma 


A.l 


indeed shows that x[^\t) >: X^^^t). 


(ii) This case follows using the exact same reasoning as for case (i) and we thus omit the 
proof. 


B 


Proof of Lemma 




The following proof is based on the proof in m Proposition 3.2]. We first describe convergence 
in distribution for a general sequence of birth-death processes and apply these results to our 
processes of interest. 

Define A(^)(-) as a continuous-time birth-death process with state space = {oi'HO I 
0 < f < oo}, where is increasing in i and limj_,.oo a^^\i) = ®)(oo). Then let 

e*'^^(a:) = argsup{a^^^(i) : a^^\i) < x}, x G [a^^^(O), a^^^(oo)). (B.l) 

ieNo 
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Denote by and the birth-death parameters associated with state 

The infinitesimal mean and variance of the process are given by 

mW(x) = + 1) - a^^'>{e‘^^\x))'^ 

- IJ^^\e^^\x))(^a^^\e‘^^\x))-a^^\e^^\x) - 1)), (B.2) 

(cr^)*'*^(x) = A^^^(e*''^^(x)) + 1) — 

+ ^i^^\e^^'>{x))(a^^'>{e‘^^\x))-a^^'>{e^^\x) - 1))^ (B.3) 

whenever x € [a(^)(0), a(^)(oo)). 

Stone’s theorem [181 Theorem 5.1] then establishes convergence in distribution of the se¬ 
quence of Markov processes to a limiting diffusion process. 

Theorem B.l. (Stone) Let A^^\0) = and A{0) = a a.s., witha^^'> —)• a. Then the following 
two conditions are sufficient for A^^^ A as elements of D[0,oo): 

(i) beeomes dense m M as s —>■ oo; 

(ii) For every compact subinterval U of M. 

lim m^^\x) = m{x), lim (a‘^)^^\x) = a‘^{x), (B.4) 

S^OO s^oo ^ ^ 


uniformly for x € U. 

We first focus on the fast system with scaled process D^\-) and state space = {(f — 
s)/^/s I i G No}. Naturally, lims_>.oo is dense in M, by which we mean that Vx G M, Ve > 
0,3s > 0 such that inf^^^(s) |x — yj < e, thus satisfying condition (i) of Theor 

that for this state space = [s -|- x^/sj , x G [—\/s, oo). 

Now, use that for the fast (and slow) system the expression a^^\e^^\x) 3- 1) — (x)) 

in (B.2) and (B.3) is equal to l/\/s to obtain the infinitesimal mean and variance of the process 

^H^as 


B.l 


Note 


(s) / \ 

[X) 

{alY"\x) 


- [s-3x\/sJ/i[j^^), X < 0, 

^(AW - s/r[j^^), X > 0, 

4(aW- 3 [s-3 X\/sJ/.iH^), X < 0, 

4(A(®)-3 s/x|j®^), X > 0. 


We use the scaling (5.2) and (^5.3^ to obtain 


(s) /' \ 
m-H (x) 


U{-/3-iV^7 + 
1 a‘[.(-/3-7). 

U,{l + i£±£va_A 

U( 2 -'^), 


s—[s+xy^J \ 

^/5 N 


X < 0, 
X > 0, 


X > 0. 


(B.5) 

(B. 6 ) 


(B.7) 

(B. 8 ) 


By requiring that lims^oo ^ (x) and lim^-^oo (o’!) (x) are finite, we indeed find that f and 7 

can be any value larger than 0. For every compact subinterval U of M, lims^oo mH^(x) = mH(x) 
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and lims_^oo (c’’!)= ^h(®) uniformly for x € U. Thus, condition (ii) of Theorem 
satished and A- 

Next we turn to the slow system with the associated scaled process Its state space 

(s) 

is equal to ^ and thus in the limit for s —>• oo also dense in M. Using the scaling as proposed 
in (5.2) and (5.3), the infinitesimal mean and variance of the process are 


B.l 


IS 


(^) / \ 
ml [X) 


I<0, 

-/3/rL, X > 0, 

A), ^<0, 

^ l (2 — -A), X > 0. 


(B.9) 

(B.IO) 


For every compact subinterval U of M, lims^oo mA{x) = mi^{x) and lirng^oo 
uniformly for x ^ U. Again, conditions (i) and (ii) of Theorem 


B.l 


= A(x) 

(s). 


are satisfied and Dl^"^ (t) 
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